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Abstract 

We present our result for the K —)• Trvr decay amplitudes for both the AJ = 1/2 and 3/2 processes 
with the improved Wilson fermion action. Expanding on the earlier works by Bernard et al. and by 
Donini et al., we show that mixings with four-fermion operators with wrong chirality are absent even 
for the Wilson fermion action for the parity odd process in both channels due to CPS symmetry. 
Therefore, after subtraction of an effect from the lower dimensional operator, a calculation of the 
decay amplitudes is possible without complications from operators with wrong chirality, as for the 
case with chirally symmetric lattice actions. As a first step to verify the possibility of calculations 
with the Wilson fermion action, we consider the decay amplitndes at an nnphysical qnark mass 
mx ~ 2m7r. Onr calculations are carried ont with Nf = 2 -|- 1 gauge configurations generated 
with the Iwasaki gauge action and nonperturbatively 0(o)-improved Wilson fermion action at 
a = 0.091 fm, m^n^ = 280 MeV, and mx = 580 MeV on a 32^ x 64 {La = 2.9fm) lattice. For the 
quark loops in the pengnin and disconnected contributions in the 1 = 0 channel, the combined 
hopping parameter expansion and truncated solver method work very well for variance reduction. 
We obtain, for the first time with a Wilson-type fermion action, that ReAo = 60(36) x 10“® GeV 
and ImAo = —67(56) x 10“^^ GeV for a matching scale q* = 1/a. The dependence on the matching 
scale q* for these values is weak. 

PACS numbers: ll.15.Ha, 12.38.Gc, 13.25.Es 
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I. INTRODUCTION 


Calculation of the K ^ mi decay amplitudes is very important to quantitatively under¬ 
stand the AJ = 1/2 rule in the decay of the neutral K meson system and to theoretically 
predict the direct CP violation parameter (e'/e) from the Standard Model. A direct lattice 
QCD calculation of the decay amplitudes for the AJ = 3/2 process has been attempted 
for a long time. Recently, the RBC-UKQCD Collaboration presented the results at the 
physical quark mass in Ref. Q, and those in the continuum limit in Ref. 0 in the physical 
kinematics, where the pions in the hnal state have hnite momenta. They used the domain 
wall fermion action which preserves chiral symmetry on the lattice. 

A direct calculation of the decay amplitudes for the AJ = 1/2 process has been unsuc¬ 
cessful for a long time, due to large statistical fluctuations from the disconnected diagrams. 
A hrst direct calculation was reported by the RBC-UKQCD Collaboration in Ref. |3| at a 
lattice spacing a = 0.114 fm and a pion mass m.,^ = 422 MeV on a 16^ x 32 lattice with 
the domain wall fermion action. They also presented a result at a smaller quark mass 
(?T7,jr = 330MeV) on a 24^ x 64 lattice with the same fermion action at Lattice 2011 0. 
In these two calculations, the kinematics was a K meson at rest decaying to two zero mo¬ 
mentum pions at an unphysical quark mass satisfying rriK ~ 2/71,^. The RBC-UKQCD 
Collaboration has since been attempting a direct calculation in the physical kinematics at 
the physical quark mass by utilizing G-parity boundary conditions. Their preliminary result 
was reported at Lattice 2014 0. 

An aim of the present article is to report on our calculation of the JJ —)■ tttt decay ampli¬ 
tudes with the improved Wilson fermion action for both the AJ = 1/2 and 3/2 processes. 
That such a calculation is feasible stems from a realization, as shown in the present article, 
that CPS symmetry 0 and its extensions 0 ensure that mixings with four-fermion opera¬ 
tors with wrong chirality are absent even for the Wilson fermion action for the parity odd 
process in both channels. A mixing to a lower dimension operator does occur, which gives 
unphysical contributions to the amplitudes on the lattice. However, it can be nonpertur- 
batively subtracted by imposing a renormalization condition 0, 0 . After the subtraction 
we can obtain the physical decay amplitudes by the renormalization factor having the same 
structure as for the chiral symmetry preserved case. A potential advantage with the Wilson 
fermion action over chirally symmetric lattice actions such as the domain wall action is that 
the computational cost is generally smaller. Hence, with the same amount of computational 
resources, a statistical improvement may be expected with the lattice calculation of the 
decay amplitudes, albeit this point has to be verihed by actual calculations. 

In the present work, we consider the decay of K meson to two zero momentum pions at 


an unphysical quark mass ttik ~ 2m.n-, as in Refs. [3|, 1^, as the hrst step of a study with the 
Wilson fermion action. Our calculations are carried out on a subset of gauge conhgurations 
previously generated by the PACS-CS Collaboration with the Iwasaki gauge action and the 
nonperturbatively 0(a)-improved Wilson fermion action for N^ = 2 + 1 havors at [3 = 1.9 


on a 32^ x 64 lattice [10[. The subset corresponds to the hopping parameters Kud = 0.13770 


for the up and the down quark and Kg = 0.13640 for the strange quark. We further generate 
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gauge configurations at the same parameters to improve the statistics. The total number 
of gauge conhgurations used in the present work is 480. The parameters determined from 
the hadron spectrum analysis are a = 0.091 fm for the lattice spacing, La = 2.91 fm for 
the lattice size, = 275.7(1.5) MeV and tuk = 579.7(1.3) MeV for the pion and the K 
masses. The energy of the two-pion state is shifted from 2m.,r by the two-pion interaction on 
the lattice. The energy difference between the initial K meson and the hnal two-pion state 
takes a nonzero value, AE = 21(3) MeV for the 1 = 2 channel, and 36(18) MeV for the 1 = 0 
channel on these conhgurations. In the present work we assume that these mismatches of 
the energy give only small effects to the decay amplitudes. 

This paper is organized as follows. The 77 —)■ tttt decay amplitudes can be calculated from 
the product of the K —)• tttt matrix elements of the AS = 1 four-fermion weak interaction 
operators and the Wilson coefficient functions for the operator product expansion. In Sec. UTI 
these four-fermion operators are introduced and the operator mixing among them for the 
Wilson fermion action is discussed. In Sec. m we describe the method of calculation used 
in the present work. The simulation parameters are also given. We present our results 
in Sec. IIVI and compare them with those by the RBC-UKQCD Collaboration and the 
experimental values. Conclusions of the present work are given in Sec. |Vl 

Our calculations have been carried out on the PACS-CS computer and T2K-Tsukuba at 
University of Tsukuba, the K computer at the RIKEN Advanced Institute for Computa¬ 
tional Science, SR16000 at University of Tokyo, and SR16000 at High energy Accelerator 
Research Org anization (KEK). Our preliminary results have been reported at Lattice 2013 
and 2014 
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II. A5 = 1 OPERATORS 


A. AS = 1 weak operators in the continuum theory 


The effective Hamiltonian of the K ^ Tin decay in the continuum theory can be written 


as 
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G 
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^ ^ ^ Qii.L) , 


( 1 ) 


2=1 


with r = — {y^gVtd) / ( 14 * 5 Kirf), and Zi(/j,) and yi{fi) (z = 1 , 2 , • • • 10 ) are the coefficient func¬ 
tions at renormalization scale fi. Here we consider the case /i < rUc, where three light 
quarks, up, down and strange, are the active quarks in the theory. The ten operators Qi^n) 
{i = 1,2, •••10) denote the AS = 1 four-fermion operators renormalized at p, which are 
given by 


Qi = {M){uu)ll , 

Q2 = (s X d){u X u)ll , 
Qs = (sd) ^{qq)LL , 

9 


( 2 ) 

( 3 ) 

( 4 ) 
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Q 4 = (s X d) ^(g X q)LL , 

q 

(5) 

Qz = (sd) ^{qq)LR , 

q 

( 6 ) 

ge = (s X d) ^(g x g)^^ , 

q 

(7) 

g? = ^ {sd)'^eq{qq)LR , 

q 

( 8 ) 

3 

Q8 = -{s X d)Y^ eq{q X q)LR , 

q 

(9) 

3 

gg = 2 i^d,)'^eq{qq)LL , 

q 

( 10 ) 

3 

gio = X d)Y^ eq{q x q)LL , 

( 11 ) 


<? 


where (sd){uu)L^R/L = (s 7 aj (1 — 75 )<^)(h 7 ;i(l ± 75 )^), and x means the contraction of color 
indices according to (s x d)L{u x d)L = J2ab(^a.^b)L{uhda)L- The summation for q is taken 
for the active quarks (g = u,d,s) and the electric charge takes values = +2/3 and 
ed = es = -1/3. 

In the four-dimensional space-time these operators are not all independent, satisfying the 
relations : 


Qi — —Qi + Q 2 + Qz 5 ( 12 ) 

Q9 = (3Qi-g3)/2, (13) 

Qio = (3^2 — = Q 2 + {Qi — Qz)/‘^ ^ (14) 

due to the Fierz identity. In general dimensions, however, these relations are not valid. 
Therefore, if we adopt the dimensional regularization for regularization, we should deal with 
all operators Q* for i = 1, 2 , • • • 10 as independent. 

B. Operator mixing for the Wilson fermion action 

The matrix elements calculated on the lattice are converted to those in the continuum by 
the renormalization factor for the operators. In this section we discuss the renormalization 
factor in the case of the Wilson fermion action. 

As already mentioned in Sec. Ill Al the 10 four-fermion operators Qi are not independent, 
and they may be arranged into 7 linearly independent combinations according to the ir¬ 
reducible representation of the flavor SU(3 )l x SU(3)r symmetry group. The 7 operators 
consist of (27, 1 ) + 4 • (8, 1 ) + 2 ■ (8, 8), whose components are given by 

(27, 1 ) g; = 3gi + 2g2 - ga, (is) 

(8, 1 ) Q '2 = 2gi — 2 Q 2 + Qz 5 
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^3 = 

Qb ) 

(8,8) Qy, 


SQi + 3Q2 + Qs 
Qg ) 

Qs ■ 


(16) 

(17) 


The operators Q'l 2 3 are the LL type four-fermion operators and Q 5 , 6 , 7,8 are of LR type. 

If the chiral symmetry is preserved in the regularization, mixings between operators in 
different representations are forbidden. For the Wilson fermion action, however, chiral sym¬ 
metry is broken to the vector subgroup, SU(3 )l x SU(3)yj — >■ SU(3)y. Hence mixings among 
different representations is in general allowed, and new operators arise through radiative 
corrections. However, we show below that such a problem is absent for the parity odd part 
of the operators in the list of (l^- flTT|) or of ffT5|) - (IT7)) for the Wilson fermion action employed 
in the present work. 

To investigate the operator mixing, we exploit the full set of unbroken symmetries for the 
Wilson fermion action, namely flavor SU(3)v^, parity P, charge conjugation C, and CPS 
which is the symmetry under CP transformation followed by the exchange of the d and the 
s quarks. All operators in the list (IT^ - ffTT)) are CPS = -|-1 operators, but the following 
operators also have the same quantum numbers including CPS, 

Qx = {M){dd — ss)sp+ps , 

Qy = {s X d){d X d — s X s)sp+ps , 

where {M){dd)sp+ps = {sd)s{dd)p+{M)p{dd)s and (sd)^ = sd and (sd)p = sj^d. Therefore 
we have to consider mixings including these operators. 

Let us discuss the problem in two steps, hrst considering mixings through diagrams in 
which gluons are exchanged between quarks of the four-fermion operators (gluon exchange 
diagrams), and second through penguin diagrams in which a pair of quarks from the four- 
fermion operators forms a quark loop. 

For the hrst type of mixings, it was shown in Ref. that the parity odd part of the 
LL and LR type operators, and the SP + PS type operator do not mix with each other. 
One can prove this through the use of CPS, CPS' and CPS" symmetries which holds for 
the gluon exchange diagrams, where S' and S" are the havor switching for a four-fermion 
operator (^iV’ 2 )(V’ 3 ^ 4 )rir 2 or {'tpi x ^/' 2)(^3 x ^/’ 4 )rir 2 dehned by 


(18) 

(19) 


S' : V’l ^ ^2 , ^3 ^ V’4 , 

S" : -01 y-)- , ^2 •H- -03 ■ 


( 20 ) 

( 21 ) 


The parity odd part of the LL and LR type operators in flT^ - flTT)) . which are of —VA — AV 
and VA — AV type, and that of Qx,y in ffT8|l - flT^ , which is of S'P-I-PS' type, are eigenvectors 
of the CPS' and CPS" symmetry with a different set of eigenvalues. 


( 22 ) 




CPS' 

CPS' 

LL\p=_i — 

-VA - AV 

+1 

+1 

LR\p=_i = 

VA-AV 

+1 

-1 


SP + PS 

-1 

-1 
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Therefore, Qx,y (the SP + PS type) do not mix with the operators (the LL and 

the LR type). 

Furthermore, the operators G (8, 8) (the LR type) do not mix with the LL type 
operators (Q123 ^ (27,1), (8,1)), or with Q ^ q G (8,1) (the LR type) because the gluon 
exchange diagrams do not change the flavor structure. 

In addition, the mixing between the (27,1) and (8,1) representations is forbidden by the 
flavor SU(3)\/ symmetry. To sum up, the matrix of the renormalization factor for the gluon 
exchange diagrams has the same structure as in the chiral symmetry preserved case. 

Next we investigate the possibility of unwanted mixings though the penguin diagrams. 
In the penguin diagrams for (^ 7,8 £ (8,8), a cancellation of the quark loop at the weak 
operator occurs between the d quark and the s quark contributions. This can be seen as 
the follow. The penguin diagram for the parity odd part of the operators Qj, except for the 
contribution from the spectator quarks, can be written as 

C 7 = CvA — Cav , (23) 

where 

C'rir 2 = T[s(X) (s(i)(-uM — — ss/ 2 )rir 2 ( 3 ^) d{Y)] , (24) 

at the space-time position x, with the external s quark s(X) at X and the d quark d{Y) at 
Y. Rewriting with the quark propagator Qq{x, y) for the quark q, we obtain, 

CriFa = +Qs{X,x)TiQd{x,Y)TY[{-Q^{x,x)+ Qd{x,x)/2 + Qs{x,x)/2)T2] 

-Qs{X, x) TiQd{x, x) T 2 Qd{x, Y)/2 - Qs{X, x) T 2 Qs{x, x) TiQd{x, Y)/2 . (25) 

Using the isospin symmetry = Qd-, C-j can be written by 


= -Qs{X, x) 7^ Qd{x, Y) Tr [iQd{x, x) - Q,(x, x)) 7^75 ] /2 

-Qs(X, x) 7/2 Qd{x, x) 7/275 Qd{x, Y )/2 - Qs{X, x) 7/275 Qs{x, x) 7/2 Qd{x, Y )/2 
- - Qs{X,x)'y^'y^Qd{x,Y)i:i[{Qd{x,x) - Qs{x,x ))/2 

-Qs{X, x) 7/275 Qd{x, x) 7,2 Qd{x, Y )/2 - Q^iX, x) Qs{x, x) 7/275 Qd{x, Y )/2 
= -Qs{X, x) 7/2 Qd{x, Y) Tr [Qds{x, x) 7/275 ] /2 
+Qs{X, x) 7/275 Qdix, Y) Tr [{Qds{x, x) 7/2 ] /2 
-Qs{X, x) jf,Qds{x, x) 7/275 Qd{x, Y )/2 

+Qs{X, x) 7/275 Qds{x, x) 7/2 Qdix, Y )/2 , ( 26 ) 


where Qdsix, x) = Qdix, x) — Qsix, x). Here we see that a cancellation of the quark loops at 
the weak operator occurs between the d quark and the s quark contributions. We can see 
this cancellation also in the operator Qg. 

Since this cancellation means that the renormalization factor coming from the penguin 
diagram is proportional to the quark mass difference (m^ — m^), mixings to four-fermion 
operators are absent due to the dimensional reason. In addition, the operator arising from 
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the penguin diagrams should have the flavor structure {sd){uu + dd + ss), which is different 
from that of Thus, operator mixings from (^ 7,8 ^ (8,8) to the other representations 

and their reverse are absent. These statements also hold for Qx,y in ffT8|) - ffT^ for the same 
reason, and the operators Qx,y are fully isolated in the theory. Further mixing between 
the (27,1) and (8,1) representations in the penguin diagrams is forbidden by the flavor 
SU(3)y symmetry. This concludes the proof on the absence of unwanted mixings among the 
parity-odd part of dimension 6 operators. 

Up to now, we have shown that the matrix of the renormalization factor for the parity 
odd part of the four-fermion operators in ffT^ - flTTD have the same structure as that in the 
chiral symmetry preserved case. Here we consider the mixing to lower dimensional operators. 
From CPS symmetry and the equation of motion of the quark, there is only one operator 
with the dimension less than 6, which is 

Qp = {md - ms) ■ s'y^d . (27) 

This operator also appears in the continuum, but does not yield a nonvanishing contribution 
to the physical decay amplitudes, since it is a total derivative operator. However, this is not 
valid for the Wilson fermion action due to chiral symmetry breaking by the Wilson term, 
and the operator fl27l) does give a nonzero unphysical contribution to the amplitudes on 
the lattice. This contribution should be subtracted nonperturbatively, because the mixing 
coefficient includes a power divergence due to the lattice cutoff growing as 1/a^. In the 
present work, we subtract it by imposing the following condition jsl, l9| , 

(0| g, |iF) = (0| Q, - A • Qp \K) = 0, (28) 

for each operator Qi in ([2])-([6]). The matrix of the renormalization factor of the subtracted 
operators has the same structure as in the chiral symmetry preserved case. 

Here, we mention an ambiguity in the subtraction procedure. Instead of strictly de¬ 
manding the subtraction condition (l28ll . we can choose a different subtracted operator, 
Qi = Qi P A ■ Qpy where A is a finite constant depending on the quark masses. The con¬ 
stants do not include the power divergence, and they vanish in the chiral limit. In general, 
such a finite ambiguity seems to remain in the final results of the decay amplitude for finite 
quark masses, as pointed out in Ref. . Our case, however, is not a such case for the follow¬ 
ing reason. The operator Qp can be written as Qp = {rrid — ms)/{rnd + ms) ■ — aX^) 

from the relation of the partially conserved axial vector current for the Wilson fermion ac¬ 
tion, where is the renormalized axial vector current and is the dimension 5 operator 
whose matrix element vanishes in the continuum limit. Thus, a A ' Qp term yields a con¬ 
tribution of form Ap ■ C — {7nT\aX a\K) ■ D to the decay amplitude with finite constants C 
and D, where Ap is the momentum difference between the initial and the final state. These 
contributions do not include any power divergent parts. Thus, by taking Ap —)• 0 and the 
continuum limit, we can safely estimate the physical value of the decay amplitudes without 
suffering from the ambiguity. 
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III. METHOD 


Simulation parameters 


Our calculations are carried out on a subset of gauge configurations previously generated 
by PACS-CS Collaboration with the Iwasaki gauge action and nonperturbatively 0(a)- 


improved Wilson fermion action at /3 = 1.9 on a 32^ x 64 lattice 10|. The subset corresponds 


to the hopping parameters Kud = 0.13770 for the up and the down quark and Kg = 0.13640 for 
the strange quark. In order to improve the statistics we further generate gauge conhgurations 
by two runs of the simulation. The hrst run uses the same algorithm as employed at the 


same parameters in Ref. 


10 


The trajectory length is r = 1/4 and the dead/alive link 


method with random parallel translation is used. The length of MD time, i.e., the number 
of trajectories multiplied by the trajectory length r, of this run is 6, 000 units as compared 


to 2,000 for the original run of Ref. [l0|. The second run does not use the dead/alive link 


method. All links are active, the trajectory length equals r = 1, and the length of run is 
also 6, 000 MD time units. We measure hadron Green’s functions and the decay amplitudes 
at every 25 MD time units for both runs. The total length of the run is 12, 000 MD time 
units and the total number of gauge configurations employed for the measurement is 480. 

We estimate statistical errors by the jackknife method with bins of 10 conhgurations (250 
MD time units). The parameters determined from the spectrum analysis are a = 0.091 fm 
for lattice spacing. La = 2.91 fm for spatial lattice size, and = 275.7(1.5) MeV and 
rriK = 579.7(1.3) MeV for the pion and the K meson masses. 

In the present work, we consider the decay in the unphysical kinematics, where the 
K meson decay to two zero momentum pions. The energy difference between the initial 
K meson and the hnal two-pion state is Ai7 = rriK — = 21(3) MeV for 1 = 2 and 

36(18) MeV for / = 0 on our conhgurations as shown in the following section. In the present 
work, we assume that these mismatches of the energy give only small ehects to the decay 
amplitudes. 


B. Time correlation function for K tttt 

We extract the matrix element {K\Qi\n7i] I) from the time correlation function for the 
77 —)■ TTTT process, 

1 - 

^iit) = 7^ + 5) Qi{t + S) IT,f^(G + 5,G + 1 + <5) |0) . (29) 

( 5=0 

Let us describe various features of this dehnition one by one. Firstly, Qi{t) is the subtracted 
weak operator at the time slice t dehned by 

Qiit) = , (30) 

X 

with the subtracted operator Q^{:K,t) at the space-time position (x, t) dehned in fl28|l . 








Secondly, the operator WKo{t) is the wall source for the meson at the time slice t, 


WK0{t) = -Wd{t)j5Ws{t) , 


( 31 ) 


with the wall source for the quark q = u,d, s, 


X 

(32) 

• 

(33) 


X 


We adopt = —d'y^s as the neutral K meson operator, so our correlation function has an 
extra minus from the usual convention. 

Thirdly, the operator Wl^{ti,t 2 ) in (|^ is the wall source for the two-pion state with the 
isospin J, which is dehned by 


WlfihM) = [{WAh)WAh) + W^+{t^)W^-{t2))/V^+{h^t2)\ / 2 , ( 34 ) 

wlfihM) = \{-wAti)WAh)/V 2 + V2W^+{h)w^-{t2))/V^+{h ^ ^ 2 )! / 2 , ( 35 ) 


where ITVi(t) is the wall source for vr* meson at the time slice t. 


= , ( 36 ) 

wAt) = {Wu{tAWA)-wAWd{t))/^^2 , ( 37 ) 

W^-{t) = Wu{tAWA) • ( 38 ) 


The wall source of each pion is separated by one lattice unit according to ti = and 
^2 = fjr + 1 in fl2^ to avoid contamination from Fierz-rearranged terms. 

We impose the periodic boundary condition in all directions. The summation over 6, 
where T = 64 denotes the temporal size of the lattice, is taken in fl2^ to improve the 
statistics. The time slice of the K meson is set at tx = 24 and that of the two-pion at 
= 0. The gauge conhgurations are hxed to the Coulomb gauge at the time slice of the 
wall source tx + d, ti + 5 and t 2 + 5 for each 5. 

In the calculation of the the mixing coefficient of the lower dimensional operator, we 
rewrite the subtraction of the lower dimensional operator fl28|) as 

Qi = Qi ~ fdi ■ Qp = Qi — Oii ■ P , (39) 

by fl27D . with P = sj^d and a* = {md — ms) ■ jSi. The mixing coefficient cxj in fl3^ is obtained 
from the following ratio of the time correlation function, 

T-l T-1 

Q!j = ^^{0| Wxo{tx + (5i) Qi{t + (5i) |0) j ^^(0| Wxo{tx + <^ 2 ) P{t + ^ 2 ) |0) , (40) 

(5i=0 ^ 2=0 

in the large tx — t region, where P{t) = T’(x, t). 
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C. Quark contractions for K —)• vrvr and if —)• 0 


In Fig. [T]we list all of the possible quark contractions for the K —)■ tttt time correlation 
function Gf (f) in fl29|) . Again there are a number of features, so let us describe them one by 
one. 

1. Time runs from right to left in the diagrams. 

2. There are four types of contractions labeled I, II, III and IV. 

3. The diagrams show the quark contractions for the four-fermion operator 


Q = r, r2 V'J) Um . ( 41 ) 

with the color indices a, b, c, d, where the spin matrix T 1^2 and the color matrix Tabcd 
are given, depending on the operator Qi, as 

Tl = 7/x(l — 75 ) ; r 2 = 7 ;i(l— 75 ) for Qi_2,3,4,9,10 5 (42) 

Tl = 7ai(1 “ 75 ) ) r2=7Ai(l+75) for (55,6,7,8 ) (43) 

Tabcd = dabbed for ( 5 l, 3 , 5 , 7,9 ; (44) 

Tabcd daddeb for (52,4,6,8,10 • (45) 


4. In the diagrams, unmarked line segments represent quark propagators for the u or the 
d quark, while those marked by “s” are for the strange quark. The filed circles stand 
for the wall sources for the K meson or pions. The open circles refer to the matrices 
Ti or r 2 . The trace for the spin is taken along closed quark lines. 

5. The subscript 1 and 2 attached to the four contraction types I though IV refers to the 
number of the trace for the spin. 

6 . The superscript “s” for the contractions IIIi ,2 and IVi ,2 means that the quark loop at 
the weak operator is for the strange quark. 

7. It should be noted that the location of Ti and r 2 for the contraction III® and IV® are 
switched from those for IIIi and IVi. 

8 . For the contraction IV* and IV® with i = 1, 2 the contribution of the vacuum diagram, 
( 0 |iF(tx)( 5 i(f)| 0 )( 0 |HV^(G,G -|- 1 ) 1 ) 0 , should be subtracted. 


Let ns write down some explicit examples. For the contraction I 2 , we have 


I 2 


[ VV Tr( IFd(x, t] ^ 2 ) 75115 (^ 2 ; X, t)F 2 )^^ 


xTr(IL5(x,f;fi)75lL5(fi;tx)75fV5(fK;x,f)Fi )^^ • Tafecd + (G •H-^ 2 ) /2 , (46) 
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with ti = and ^2 = + 1, where the trace is taken for the spin. The three types of Wq 

{q = d, s) in are the wall source propagators for the quark q dehned by 

^(47) 

(48) 


ir,(x, (; (o) = J]] Q,{x, (; y, «o) , 

y 

Wq{to; X, f) = ^ Qq{y, to; x, t) = 75 W^q(x, t; fo)^75 , 
y 

= '^Wq{^,t;to) , 

X 

with the quark propagator Qq{yi,t;y,to). Similarly, the contraction Ii is given by 


(49) 


I. = 


EE Tr J (iyrf(x, t; t 2 )-i 5 Wdit 2 ; X, t)V 2 )^,^ 

1 h r d X 

{Wd{^,f,tlh5Wd{ti;tKh5Ws{tK]^,t)'^l)da ' ^abcd + (^1 ^ ^ 2 ) 


*- a,b,c,d X 

X 


/2 ,(50) 


where the trace is taken for the spin. 

The contraction II 2 is given by 

■■'-[EE Tr( lTrf(x, t; t2)75fkrf(^2; ti)75fkrf(^i; x, t)T2)dc 

a, 6 ,c,d X 

XTT{Wd{yi,f,tKh5Ws{tK;^,t)'^l)ba''^abcd + {h ^ t2) /2, (51) 

and the contraction IIi is given by 

III = [ ^ X^Tr {Wdix,t;t2)'y5Wd{t2;ti)'y5Wd{ti;x,t)T2)hc 

a.b.c.d X 


x{Wdixd;tKh5Ws{tK;x,t)Ti)a^ 
The contraction III 2 is given by 


• Tabcd + (^1 "t-t ^ 2 ) /2 . (52) 


IIIo = 


EE Tr( Wd{x, t; t2)75i^d(t2; ^i)75i^d(^i; ^x)75i^s(^i^; x, t)ri )^^ 

*- a,b,c,d X 

X 


Tr( Qd{x, t; X, f)r2 ■ T^bcd + {h ^ ^ 2 ) 


/ 2 , 


(53) 


where the quark loop for the d quark Qd{x, t; x, t) is calculated by the the stochastic method, 
whose detail is discussed in the next section. The contraction III 2 is obtained by changing 
Qd(x., t; X, t) to the quark loop for the s quark Qsix, t; x, t). 

Having constructed various quark contractions, we can build the K —)■ tttt time correlation 
function Gf {t) for the operators Qi in the isospin channel I in the following way. For the 
1 = 2 case, we have 


G[=^ 

G\=^ 

^1=2 

^8 



(I 2 

(I 2 

(I 2 


Ii) = G 


1=2 

2 


Ii) 


1 


Ii) 





2 

3 ’ 


(54) 

(55) 

(56) 
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where ri ^2 and Tabcd in each contractions should be chosen according to fH2|) - fH5l) for each 
operator. The relation among different operators (I5T|) follows from the Fierz identity. 

The formulae for the 1 = 0 channel are given as follows: 

for {sd){uu) = Qi ^2 

G'=° = y^(-l2-2-Ii + 3-Il2 + 3-T2) , (57) 

for {sd){uu + dd + ss) = < 53 , 4 , 5,6 

= y|(-l2 + 2-Il2-IIi + 2-T2-Ti + T^-T^) , (58) 

for (sd)(-U'U — — ss/2) = <57,8,9,10 

= y| (-I 2 - Ii + II 2 + III + T 2 + Ti - T^ + Tl) , (59) 

with Tj = Illj — IVj and T^ = III| — IV® {i = 1, 2), where ri ,2 and Tabcd in each contractions 
should be chosen according to (H2D - flT5l) for each operator. 

We now turn to the quark contractions needed to subtract the contribution of the lower 
dimension operator Qp. In Fig. [2] we list all of the possible quark contractions for the K-to- 
vacuum time correlation function Gk^o = {0|lFxo(t_ft')<5i(^)|0) in fHOj) . The notations are 
the same as for Fig. [H The contractions for the operators Qi are given by 

for {sd){uu) = < 5 i ,2 

Gk^o = — V 2 , (60) 

for {sd){uu + dd + ss) = < 53 , 4 , 5,6 

Gk^o = -2 ■ V 2 + Vi - V® + (61) 

for {sd){uu — dd/2 — ss/2) = < 57 , 8 , 9,10 

Gk^, = (-V 2 - Vi + V® - V®)/2 , (62) 

where Fi ,2 and Tabcd in each contractions should be chosen according to (H2ll - (H5ll for each 
operator. We can obtain the mixing coefficient of the lower dimensional operator o;* in (l39|) 
by dividing these by the time correlation function {0\WKoitK)Pit)\0) as (j40]). 

The K 7171 time correlation function for the operator P = sj^d is calculated by 

T—1 

GUt) = ^ 5^(0| WK0{tK + S) P{t + 5) + + 1 + |0) = -^ Tp , (63) 

where Tp = IIIp — IVp, and the contractions IIIp and IVp are shown in Fig. [3l For the 
contractions IVp, the contribution of the vacuum diagram {0|iF(tp:)P(t)|0)(0|hF/,,.(f^,+ 
1)|0) is supposed to be subtracted. 

The K ^ 7171 time correlation for the subtracted operator <5^ = <5i — cii • -P is given 
by subtracting the contributions of Oj ■ P from those of the operator Qi. We write these 
subtractions as 

III ^ III - tti ^ IIIp , (64) 

\/6 

IV ^ IV - cii ^ IVp , (65) 

\/6 
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dividing into the connected (III and IIIp) and the disconnected contractions (IV and IVp), 
where III means the total contribution from the two contractions Illj and Illf for i = 1,2, 
and, similarly, for IV. 


D. Calculation of quark loop 


The quark loop at the weak operator Q{x,x), i.e., the quark propagator starting from 
the position of the weak operator and ending at the same position, appears in the quark 
contractions III, IV, IIT and IV^ for the it' —)■ tttt process, and V and V^ for the K ^ 0 
process, as shown in the previous subsection. We calculate them by the stochastic method 
with the hopping parameter expansion technique (HPE) and the truncated solver method 


(TSM) proposed in Ref. [13 


The Wilson fermion action can be written as 


= = - D)^/J = ^|JM{1- D)^/J , (66) 

where D = M~^D and 

{M‘il;){x) = (1 - kCsw{o- ■ F{x))/2) ‘ipix) , (67) 

(Di,)(x) = (£>J + D;) V-) (x) , (68) 

{D+ij){x) = k{1- 7^) Uf,{x)'ip{x + fi) , (69) 

{D-i)){x) = K (1 + 7^) Ul{x - ^)i){x - fx) . (70) 


The quark held of the Wilson fermion ^|J is related to that in the continuum theory by 
= \/^ • xf) in the tree order. From fl66l) the quark propagator Q can be written by a 
hopping parameter expansion form as 

oo k—1 

Q = ^ ^ (71) 

n=0 n=0 

for any integer value of k. We use this expansion to calculate the quark loop Q{x,x) at the 
weak operator. In this case, the terms with odd power of D do not contribute, thus 


Q{x, x) = {M-^ + + D^W-^){x, x) , 

for k = A. We can replace term by 

= 1 ] {d*d; + d;d*) , 


(72) 


(73) 


with = M Using these expressions, we calculate the quark loop by the stochastic 

method according to 


Nr 


Q{x,t;x,t) = • 


(74) 


2=1 
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The function S'i(x, t) is defined by 


^i(x,f) = ^(M ^ + D2lM ^ + D^W ^)(x,f;y,f) , 

y 

where we introduce an U{1) noise ^j(x, t) which satisfies 


^^(x-y) 


1 


Nr 

5^C(x,t)6(y,^) 


i=l 




(75) 


(76) 


for Nr —)■ cx). The effect of HPE for the quark loop is to remove the D and terms in 
fl75p explicitly which make only statistical noise. We hnd that HPE reduces the statistical 
error of the decay amplitudes to about 50% compared with the normal stochastic method. 
We also implement the truncated solver method (TSM) for fl7^ by 


Nr+Nt 

^ , (77) 

i=NR+l 

where Sf (x, t) is a value given with the quark propagator W~^ in (175|l calculated with a 
loose stopping condition, and S'j(x, t) is that with a stringent condition. We set A% = 5 and 
the stopping condition R = |hPa; — .^|/|.^| < 1.2 x 10“® with x denoting the iterative solution 
of Wx = ^ for S'f(x,f), and = 1 and R < 10 for S'i(x, f) in the present work. The 
numerical cost of TSM (iTTll is about twice of that without TSM (ff4ll with Nr = 1. 


Nr 


Q(x, X, t) = — V C(x, t) [^i(x, t) - Sj (x, t)] 


i=l 


E. Time correlation function for Trvr —)• tttt 

We calculate two types of time correlation functions for ttyt —>■ ttyt to obtain the nor¬ 
malization factors which are needed to extract the matrix elements (iPlQjTrvr;/) from the 
time correlation function Gf(t) in fl29|) . These are point-wall and wall-wall time correlation 
functions, which are dehned by 

1 

Gpw{t) = 7^ X^(0| (5, 1 h) |0) , (78) 

5=0 

1 

Gwwi^) = j + 5, t + 1 + 5) -F (5, 1 -h 5) |0) , (79) 

5=0 

where (7r7r)'^(t) is the operator for the two-pion system with the isospin /, 

(7r7r)^=2(t) = ^(^vr°(x,t)7r°(y, t) + 7r+(x,t)7r"(y,t))/\/3 , (80) 

x,y 

(7r7r)^=°(t) = ^(-vr°(x,t)7r°(y, t )/\/2 + \/27r+(x,t)7r"(y, t)^/V3 , (81) 

x,y 
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with the operator for tt* meson 7r*(x, t) dehned by 


7r+(x,t) = -(J(x,t) 75 M(x,t) , (82) 

7r°(x, t) = (h(x, t) 75 M(x, t) - (i(x, t) 75 (i(x, t)) /\/2 , (83) 

TT" (x, t) = m(x, t) 75 d(x, t) . (84) 


The operator hl/)(^(ti,^ 2 ) (-f = 0,2) is defined by (l34)) and (l35D . In the present work we set 
fjr = 0 in fITSli and fl79l) . 

In Fig. m we list all of the possible quark contractions for the time correlation function for 
the TTvr TTTT processes. Time runs from right to left in the diagrams. There are four types 
of contractions, D, C, R and V. The hlled circles represent the wall source in Il36|) - (l38|) 
or the point source vr* in fl82l) - fl84ll for the pion. For the contraction V the contribution of the 
vacuum diagram, {0|(7r7r)'^(f)|0){0|IF^^(f^, f^ + l)|0) or (0|IT,(^(f, f+l)|0)(0|IF)(^(f^, t^ + l)|0), 
is supposed to be subtracted. 

For example, the explicit form for the contraction C for the point-wall time correlation 
function Gp-^{t) in fl78|) is given by 


C 


[ ^ Tr( Wditi^y, t) 75 W^<i(y, i] t 2 )lbWd{t 2 ] x, f) 75 lFrf(x, f; ^ 1)75 ) 
x,y 


-|- {ti O ^2) 


/ 2 , 


(85) 


with ti = and ^2 = tn + ^, and that for the wall-wall correlation function G(^/^(t) in fl79|l 
by 


C 


Tr( Wditi; U)^5Wdit4; t2)75Rd(^2; i3)75Rd(^3; ^ 1)75 ) 


+ {ti -f-)- ^ 2 ) + (^3 77 ^ 4 ) + (ti 77 t2, ts 77 ^ 4 ) 


/4, ( 86 ) 


with ti = t 2 = + 1, ts = t, = t + 1 , where the trace is taken for the color and the 

spin indices. The wall source propagators Wd are dehned by (07])-(09]). 

For the calculation of the contraction R for the point-wall time correlation function 
Gpw{t) in (178]), we use the stochastic method according to 


R = 


Nr 


Tr( Wd{ti, t2)'y5Wd{t2-, X, f )75 Z,(x, t)C(y, t) 75 lFd(y, t; ^ 4)75 

R ,•_1 ,, 


2=1 x,y 


+ {ti 77 ^2) /2 , 


(87) 


with ti = and ^2 = +1, where ^i(y, t) is an 7/(1) noise which satishes fl76]) . and .^i(x, t) 

is dehned by 

Zi(x, t) = Y^ hF"^(x, t; y, t) Ci(y, t) , (88) 

y 

with the kernel of the Wilson fermion W in (166]) . The contraction V is also calculated 
by using Zj(x, f). In actual calculations, we hnd that relaxing the stopping condition to 
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\Wx — ^l/l^l < 1.2 X 10“® for the calculation of makes only negligible effects to 

the hnal result, compared with the statistical error. Thus, we adopt this loose stopping 
condition with = 6 in flHTl) . 

The quark contraction for the time correlation function, Gp^^it) or given by 

= D - C , (89) 

G'^=o = D + ic-3R+^V . (90) 


IV. RESULTS 

A. Time correlation function for Trvr —)■ ttyt 

Fig. O shows the contributions of the four types of contractions, D, C, R, V, for the time 
correlation function for tttt —)■ tttt, with those for the point-wall function G^p^/it) in fl78|l 
plotted on the left and those for the wall-wall function Gly^^it) in fl79|) on the right panel. 
The source operator is placed at = 0. The time correlation functions behave in the large 
time region as 

Gpw{t) = -t- j -h G^ , (91) 

Gwwit) = Kn ■ , (92) 

where is the energy of the two-pion system with the isospin /, is a constant whose 
form is irrelevant, and 

Ai, = (0|fFi(0,l)|7r7r; J)V(7r7r;/|7r7r;/) . (93) 

The constant terms G^ and in fl9Tl) and fl92|) come from the two pions propagating in the 
opposite time directions (he., around-the-world effect for the two-pion operator). 

The effective mass of the point-wall time correlation function Gpy^it) is plotted in Fig. El 
where the effective mass fn^ff at t is given by 

Gpwjt + 1) - Gpy^rjt + 4) ^ f{t + l] rrieff) - /(t -h 4; nieff) 

Gpwit) - Gpwit + 3) fit] rueff) - f{t + 3; nieff) 

with the function /(t; nieff) = exp(—rUe// • t) + exp(—me// • (T — t)). We hnd plateaus in 
the time region t > 9 for both / = 0 and 2, albeit admittedly much noisier for / = 0 than 
for 1 = 2. Compared with the value 2m,r plotted in blue, the two-pion energy for 1 = 2 
is larger, signifying repulsive interaction of the two-pion system, whereas that for J = 0 is 
smaller showing attractive interaction. 

In the extraction of the matrix elements (iF|Qj|7r7r; I) from the time correlation function 
Gj (t) in ([29]), the values of and are needed. Since the statistical error of the point- 
wall correlation function Gpyy{t) is smaller than that for the wall-wall function G^-^^{t), we 
hrst extract the energy from Gp^^{t), and then extract the amplitude A^^ from G\y^{t) 
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by fitting to fl9^ with the determined value of and regarding and as unknown 
parameters. The results for and are 

Eif = 0.2567(14) , Aif = 2.513(27) x lO^o , 

Ei=^ = 0.2499(83) , = 2.41(13) x 10^° , (95) 

in the lattice unit, where we adopt the htting range t = [9, 32] for 1 = 2 and t = [9,12] for 
J = 0. 

The mass of the pion and the K meson obtained in the present work are = 0.12671(71) 
and irtK = 0.26641(58) in the lattice unit. The energy difference between the initial K meson 
and the hnal two-pion state, = rriK — E^^, is AE^^^^ = 0.0097(14) (21(3) MeV) and 
AE^^^ = 0.0165(83) (36(18) MeV). In the present work, we assume that these violations of 
energy conservation yield only small effects to the results for the 77 —)■ tttt decay amplitudes. 


B. Time correlation function for K —> vrTT in the 7 = 0 channel 

In Fig. [7] we demonstrate the effects of the truncated solver method. The four panels (a)~ 
(d) show the contributions of the contractions III and IV to the time correlation functions 
for Q 2 and Qg at f = 9. In each panel, the data at x = 0 shows the result of a stochastic 
estimate with a stringent stopping condition, while those at x = 1, • • •, 6 are obtained with 
a loose stopping condition, with an identical noise vector employed for x = 0 and x = 1. 
Thus, the difference between the data at x = 0 and x = 1 corresponds to the hrst term of 
(1771) for Nr = 1, and the data at x = 2, • • •, 6 (AV = 5, Nr + Nr = 6 ) correspond to the 
components of the second term of fl77|l . We hnd that the hrst term is negligible compared 
with the statistical error for all channels. Thus we can neglect it, and estimate the quark 
loop contribution by only the second term as 

Nt+Nr 

Q(x,t;x,t) = ^;(x,t) ^f(x,t) . (96) 

1=1 

The contribution given by the sum fl96l) is plotted at x = 7 in each panel. We see that the 
statistics is signihcantly improved by increasing the number of random numbers from 1 to 

6 . 

The results for the 7 = 0 77 —)■ tttt time correlation function for the operator Q 2 ( ^^^^(t) 
in (12^ ) are plotted in Fig. [HI The time slice of the two-pion is set at = 0 and the 77 
meson at = 24, while the operator Qi{t) runs over the whole time extent as explained 
before. In the panels (a) and (b), we observe a large cancellation between the contributions 
from the operator Q 2 and the subtraction term 02 • P for both contractions III and IV. 
In panel (c) we hnd that the contribution from the contraction IV is similar in magnitude 
to that from the contraction I. This appears different from the previous work by RBC- 
UKQCD Collaboration with the domain wall fermion action in Refs, jl, |^. In panel (d), 
we compare the correlation functions calculated with TSM and without TSM. We hnd that 
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TSM significantly improves the statistics. The numerical cost of TSM is about twice of that 
without TSM. Thus TSM is a very efficient method. 

The results for Qq in the 1 = 0 channel are plotted in Fig. [9l Here, also, we find a large 
cancellation between the contributions of Qq and the subtraction ■ P for both contractions 
III and IV (see panels (a) and (b)). In panel (c) a large cancellation is observed between 
the contraction I and II, which is not the case for Q 2 - An efficiency of TSM is observed also 
for Qq in panel (d). 


C. K ^ TTvr matrix elements 


In order to extract the K ^ mi matrix element, we consider an effective matrix element 
which behaves as M-{t) = M- = {K\ Qj(0, 0) Itttt; I) in the time region tx ^ t ^ t-,,. 
It can be constructed from the time correlation function Gf (t) in (12^ by 

Ml{t) = Gl{t)I^AKA{^ ■ ■ Q^K{tK-t)+EUt-tA X (-1) . (97) 


Here, the K meson mass rriK and the energy of the two-pion state are fixed at the 
values obtained from the correlation function of the K meson and the tttt —)■ mi. The 
factor (—1) comes from the convention of the operator in fl3T]) . The constant = 
{^\Wk\K)’^/{K\K) is estimated from the wall-wall propagator of the K meson, with the 
value Ak = 8.949(34) x 10® in the lattice unit. The constant is defined by (|9^ and 
its value is given by fl95|) . The dimensionless constant is the Lellouch-Liischer factor 14 
given by 




{K\K) ■ (tttt; JI tttt;/)/ 


(98) 


where V is the lattice volume V = L^, is the two-pion scattering phase shift for the 

two-pion system with the isospin I at the scattering momentum = (ii^^^)^/4 — m^, and 
(j){q) is the function defined by 


tan 0(g) = -7i^/‘^q/Zoo{l] q) , 


(99) 


with the spherical zeta function. 


Zoo{s;q) 



V 

neZ3 


( 100 ) 


at g = p{2ti/L). In the noninteracting two-pion case, the factor takes the form {F^Y = 
(Fifree)' = (2m^v) • (2m^VY/V\ 

For the 1 = 0 channel the statistics in the present work is not sufficient to obtain the 
scattering phase shift. We therefore use the factor for the noninteracting case, leaving a 
precise estimation of the factor to study in the future. For the 1 = 2 case the phase shift 
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is obtained with a sufficient statistics at the needed momentum. Because the scattering 
momentum p takes a small value, p = 2.053(97) x 10“^ (44.7(2.1) MeV) in our case, the 
phase shift can be approximated by = p{d5^^‘^{p)/dp) + O(p^). We neglect the cubic 

term, and hnd F-^"^/F|free = 0.9254(62). 

Our results for the effective matrix elements for several representative channels are shown 
in Fig. [To] and Fig.[TT]for 1 = 2, and in Fig.[T2]for / = 0, where the matrix elements calculated 
with tx = 22, 24, and 26 are plotted. We hnd plateaux for the effective matrix elements 
over the time interval t = [9, 12] which are independent of the value of tx- This means that 
the around-the-world effect of the two-pion operator is negligible in this time region. 

We extract the matrix element M/ = {K\ Qi{0, 0) Itttt; I) by a constant £t of the effective 
amplitude for tx = 24 in the time interval t = [9,12]. Our results for the 1 = 2 channel (the 
AJ = 3/2 process) are tabulated in the second column in Table HI where the relation among 
the matrix elements fl5T|) is used. The results for the 1 = 0 channel (the AJ = 1/2 process) 
are tabulated in the second column in Table HU Here we do not use the operator relations 
dia-dll, and treat each of 10 operators as independent. 


D. K —)• TTTr decay amplitudes 


The renormalized matrix elements (/i) are obtained from the bare matrix elements 
on the lattice extracted in the previous section by multiplying with the renormalization 
factors as 


10 


( 101 ) 


M. (?•) = !] A#/Xih'") ■ 

i=i 

The renormalization factors Zij{q*a) for our choice of the fermion and gluon actions have 


been calculated by perturbation theory in one-loop order in Ref. 15|. A nonperturbatively 


determination is not yet available. For the renormalization in the continuum theory, we 
adopt the modihed minimal subtraction scheme (MS) with naive dimensional regularization 
scheme (NDR). We choose two values q* = 1/a and vr/a as the matching scale from the lattice 
to the continuum theory in order to estimate the systematic error coming from higher orders 
of perturbation theory. Large tadpole contributions in the renormalization factors for the 
lattice perturbation theory are subtracted by the mean-held improvement. We use a mean- 
held improved value in the MS scheme for the coupling constant, which is given from the 
bare coupling constant by 

= {CoP + 8CiR)/g^-0.1006 + OmU9-Nf + {ll-2Nf/3)/{8Tr^)-\og{q*a) (102) 

for our gluon and fermion actions, where Co = 1 — 8Ci and Ci = —0.331 are the parameters 
in the gluon action, and P is the expectation value of the plaquette and R is that of the 
1x2 Wilson loop. The detail of the procedure was discussed in Ref. [l^. From the values 
P = 0.572059(31) and R = 0.338902(47) given in Ref. 10], we obtain g^ = 2.699 at 
q* = 1/a and g^ = 1.996 at g* = vr/a. 
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The decay amplitudes Ai (/ = 0,2) are calculated from ([T]) as 


10 


A, = 




where 


*J=i 


G 


Ci(fi) = + mW) 


( 103 ) 


(104) 


The explicit form of the functions Zi{^) and yi{^) in the NDR scheme have been given 


in Ref. 12j. The functions Uij{q*,y,) are the running factor of the operators Qi from the 


scale q* to p for the number of active quark flavors equal to Nf = 3, which have also been 
given in Ref. [l^. In the present work we set /i = rric = 1.3 GeV in fllU3p and evaluate the 
two functions Zi{fi) and i/i(/i) with the Standard Model parameters tabulated in Table HTTl 
We adopt the standard representation of the CKM matrix, in which CP violation enters 
entirely through the complex phase of thus r = — (y*^Vtd) / (y*gVud)- The values of the 
two functions are tabulated in Table HVl 

From fllOll) and fll03p the decay amplitudes can be written in terms of the bare matrix 
element M- as 


10 


10 


/1, = E = 


(105) 


2=1 


2 = 1 


where 


10 


Gj = ^ Zij{q*a)Ujk{q*,fi)Ck{fi) . 

j,k=l 


(106) 


The constant Ci should be independent of /i and g*, and depend only on the lattice cutoff 
1/a, if we work in the full order of perturbation theory. We define T* and by fll04p for 
Ci- The values of these quantities for q* = 1/a and n/a at fi = rric = 1.3 GeV are given in 
Table EVl 

Our hnal results of the decay amplitudes are given in Table |Vl The direct CP violation 
parameter e'/e is obtained by 


Re(eYe) 


u / Imyl2 ImRo \ 

\/2|e| \Re742 ReRo / 


(107) 


with uj = ReR2/Re44o, where the experimental value of the indirect CP violation parameter 
|e| = 2.228 X 10“^ is used in the estimation. The statistical errors are estimated by the 
jackknife procedure with a bin size of 10 conhgurations (250 MD time units). We also list 
results of the RBG-UKQGD Gollaboration at the similar quark masses ( tUtt = 422 MeV 
and 330 MeV ^ ) for comparison. These two cases are calculated with the unphysical 
kinematics at mx ~ 2m.n-, as in our calculation. In the table, the results of the RBG- 
UKQGD Gollaboration for the AI = 3/2 process obtained at the physical quark mass with 
the physical kinematics, where the pions in the hnal state have hnite momenta, in the 
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continuum limit presented in Ref. [2|], and the experimental values are also tabulated. We 
note that our results with an unphysical kinematics can not be directly compared with these 
values at the physical quark mass. 

From Table IVl we learn that the dependence on q* is negligible for most of the decay ampli¬ 
tudes, but it is very large for lmA 2 . A nonperturbative determination of the renormalization 
factor is necessary to obtain a reliable result for this value. We hnd a large enhancement 
of the AJ = 1/2 process over that for the AJ = 3/2 at our quark mass = 280 MeV. 
The RBC-UKQCD Collaboration found that the enhancement was explained by the follow¬ 
ing numerical mechanism [^: A large cancellation between two dominant quark diagrams 
occurs for the AJ = 3/2 process, rendering ReA 2 small, while such a cancellation does take 
place for the AJ = 1/2 process. We conhrm this numerical mechanism also in our case. 

Our result for Aq, particularly for the imaginary part, still has a large statistical error 
so that we do not obtain a nonzero result for Re(eYe) over the error. We observe that the 
results for Aq by the RBC-UKQCD Collaboration at a similar pion mass m-n = 330 MeV 
have smaller errors than ours. This is because they use a different two-pion operator for 
which the wall sources for the two pions are separated by 5 = 4 in the time direction, 
and they set the htting range closer to the two-pion source than our case in extracting the 
matrix elements from the time correlation functions. Improving statistics by devising a more 
efficient operator for the two-pion state is an important work reserved for the future. 

The contributions of the bare matrix element M/ to the decay amplitude Aj {Aj{i) in 
(11051) ) are tabulated in Table [T] for the AJ = 3/2 and in Table UTI for the AJ =1/2 process. 
We hnd that the main contribution to ReA 2 comes from the operator Qi and Q 2 , and that 
to ImA 2 from Qg. The main contribution to ReAg comes from the operator Q 2 and that to 
ImAo from Qq. 


V. CONCLUSIONS 

In the present work we have shown that mixings with four-fermion operators with wrong 
chirality are absent even for the Wilson fermion action for the parity odd process due to CPS 
symmetry. Therefore, after subtraction of an effect from the lower dimensional operator, 
a calculation of the decay amplitudes is possible without additional calculations for the 
operator with wrong chirality. This is the same situation for chirally symmetric lattice 
actions such as the domain wall action. A potential advantage with the Wilson fermion 
action over chirally symmetric lattice actions is that the computational cost is generally 
smaller. Hence, with the same amount of computational resources, a statistical improvement 
may be expected. 

As the hrst step of a study to verify the possibility of calculations, we considered the K 
meson decay amplitude for both the AJ = 1/2 and 3/2 channels with the Wilson fermion 
action at an unphysical quark mass rriK ~ 2m^. We have found that the stochastic method 
with the hopping parameter expansion technique and the truncated solver method are very 
efficient for variance reduction, yielding a hrst result for the J = 0 amplitude with the Wilson 
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fermion action. 

We have been able to show a large enhancement of the AJ = 1/2 process over that 
for the AJ = 3/2 at our quark mass (m^ = 275.7(1.5) MeV and rriK = 579.7(1.3) MeV). 
However, our result for Hq, particularly for the imaginary part, still has a large statistical 
error so that we have not obtained a nonzero result for Re(eYe) over the error. For the 1 = 0 
two-pion system, the statistics in the present work are not sufficient to obtain the scattering 
phase shift. We therefore used the Lellouch-Liischer factor for the noninteracting case in 
the calculation of the AJ = 1/2 process. Improving statistics by devising a more efficient 
operator for the J = 0 two-pion state is an important work reserved for the future. 

Our calculation is carried out away from the physical quark masses, and the decay of the 
K meson to two zero momentum pions at ttik ~ is considered. Clearly, we need to 
work toward smaller quark masses and a more realistic kinematics in which the two pions 
carry hnite momenta. This will be a major challenge that we now have to face. 
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Figures 



FIG. 1: Quark contractions for the time correlation function for the K —)■ tttt process for the 
operator Qi (i = 1, 2, • • • 10). 
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FIG. 2; Quark contractions for the time correlation function for the ii' —)■ 0 process for the 
operator Qi (i = 1, 2, • • • 10). 



FIG. 3: Quark contractions for the time correlation function for the K —)■ tttt process for the 
operator P = sj^d. 



FIG. 4: Quark contractions for the time correlation function for tttt —)■ tttt. 
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FIG. 5: Four types of contractions for the time correlation function for tttt —)• tttt. Left panel 
shows those for the point-wall function Gpyy{t) in (1751) and right for the wall-wall function 
in (17^ . 



FIG. 6: Effective mass of the time correlation function Gpyy{t) for vr-Tr — )• tttt with the isospin 
1 = 0 and 1 = 2. Twice of the effective mass for the pion is also plotted for a comparison. 


26 

















































FIG. 7: Effect of the truncated solver method. Panels (a)-(d) show the contributions of the 
contraction III and IV to the time correlation functions for Q 2 and Qq at t = 9. In each panel, the 
data at X = 0 show the contribution obtained by the usual stochastic method (|74p with Nji = 1, i.e., 
the contribution given by setting the quark loop = ^*{x,t) Si{x,t) for i = 1. The data 

at X = 1, 2, • • •, 6 correspond to the contributions given by setting Q{x, t; x, t) = 4*(^) 
for X = i = 1, 2, • • • , 6 {Nr + Nt = 1 + 5) with Sf (x, t) in ([771l . The data at x = 7 are average of 
the data at x = 1, 2, • • •, 6. 
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FIG. 8: Time correlation function for the operator Q 2 for the A/ = 1/2 K ^ tttt process, 
in (|29p . The time slices of the two-pion and the K meson are set at G = 0 and tx = 24, while the 
operator Qi runs over the whole time extent, (a) Contributions of the contraction III for Q 2 , 02 - P 
and Q 2 = Q 2 — Ci 2 • P- (b) Contributions of the contraction IV for Q 2 , 02 • T and Q 2 = Q 2 — 0 ( 2 - P- 
(c) Contributions from each type of contractions for Q 2 , (d) Total correlation functions calculated 
with TSM and without TSM. 
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FIG. 9: Time correlation function for the operator Qq for the A/ = 1/2 K ^ tttt process, Gq ^{t) 
in (j29l) . following the same convention as in Fig. [5J 
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FIG. 10: Effective matrix element of in (|97p . in the lattice unit. The time slice of the 

two-pion is set at tjr = 0. The operator Qi{t) runs over the whole time extent. The matrix elements 
given with the K meson at time slice tx = 22, 24, 26 are plotted. 



FIG. 11: Effective matrix elements of Mj i 



(I97p . following the same convention as in Fig. 1101 



FIG. 12: Effective matrix elements of in (j97|] . following the same convention as in Fig. 1101 
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Tables 


TABLE I: Decay amplitude for the AI = 3/2 process. The second column gives the bare matrix 
elements for Qi in the lattice unit. The other columns are their contribution to A 2 ( ^ 2 ( 1 ) 

in (11051) ) ioi q* = 1/a and n/a. 

q* = 1/a Q* = '^/cL 

i _M/_ ReAa (GeV) ImeAa (GeV) _ ReAa (GeV) ImeAa (GeV) 

1 2.256(35) X 10“3 -1.887(29) x 10"°® 0 -1.452(23) x 10“°® 0 

2 = M/=2 4.330(68) x 10"°® 0 3.920(61) x 10“°® 0 

7 9.85(11) X 10"2 1.053(12) x 10"^° 2.772(32) x 10"^® 3.172(36) x 10"®° 2.100(24) x 10"®® 

8 3.242(37) x 10"® -2.722(31) x 10"®° -1.670(19) x 10"®^ -4.124(47) x 10“®° -1.156(13) x 10“®^ 

9 = 3/2 • M/=2 -1.140(18) X 10"®2 3.762(59) x 10"®® 3.739(58) x 10“®^ 3.409(53) x 10“®® 
10 = 3/2 • M/=2 3.771(59) x 10"®° -1.756(27) x 10"®® 4.372(68) x 10"®° -1.409(22) x 10“®® 

Total - 2.426(38) x 10"°® -1.192(14) x 10“®^ 2.460(38) x 10"°® -7.457(83) x 10“®® 


TABLE II: Decay amplitude for the AI =1/2 process. The second column gives the bare matrix 
elements M/^° for Qi in the lattice unit. The other columns are their contribution to Aq ( Ao(i) 
in (11051) ) for O'* = 1/a and ir/a. 




q* = 1/a 

q* = vr/a 

i 

M/ 

ReAo (GeV) ImeAo (GeV) 

ReAo (GeV) ImeAo (GeV) 


1 u.oi^i.oj X 10 ^ 

2 3.6(1.4) X 10"2 

3 7.2(3.7) X 10"2 

4 1.06(40) X 10"® 

5 -1.0(4.3) X 10"2 

6 -2.0(1.1) X 10"® 

7 2.42(18) X 10"® 

8 7.46(54) X 10"® 

9 -3.0(1.4) X 10"2 

10 0.0(1.2) X 10"2 

Total 


6.8(2.8) X 10"°®’ 
-1.25(65) X 10"°® 
5.3(2.0) X 10"°® 
1.5(5.9) X 10"°° 
-8.4(4.6) X 10"°® 
2.58(19) X 10"®° 
-6.26(45) X 10"®° 
1.02(48) X 10"®® 
0.0(1.4) X 10"®® 
6.0(3.6) X 10"°®' 


0 

0 

-2.5(1.3) X 10"®® 
6.6(2.5) X 10"®® 
1.7(6.8) X 10"®2 
-1.03(56) X 10"®° 
6.81(50) X 10"®® 
-3.84(28) X 10"®2 
-3.4(1.6) X 10"®2 
-0.1(6.4) X 10"®® 
-6.7(5.6) X 10"®® 


6.2(2.5) X 10“°^ 
-1.7(8.7) X 10“°® 
6.2(2.4) X 10“°® 
1.9(7.4) X 10“°° 
-7.7(4.2) X 10“°® 
7.79(57) X 10"®° 
-9.48(68) X 10“®° 
-3.4(1.6) X 10“®® 
0.0(1.6) X 10"®® 
5.6(3.2) X 10"°®’ 


0 

0 

-2.1(1.1) X 10"®® 
6.1(2.3) X 10"®® 
1.8(7.1) X 10"®2 
-8.8(4.8) X 10"®® 
5.16(38) X 10"®® 
-2.66(19) X 10"®2 
-3.1(1.4) X 10"®2 
-0.1(5.2) X 10"®® 
-5.2(4.8) X 10“®® 
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TABLE III: Standard model parameters used to evaluate the decay amplitudes in the present 
work (from Ref. [l^). r = — {V^*Vtd) / {Vus^ud) and A^ is the lambda QCD for Nf = 5 theory. 


The standard representation of the CKM matrix of Ref. 
enters entirely through a complex phase of Vtd, thus r. 


is adopted, where the CP violation 


mz 

91.188 GeV 

mw 

80.385 GeV 

mt 

173 GeV 

mb 

4.2 GeV 

me 

1.3 GeV 

MS 

0.23135 GeV 

a ( at /X = m\Y ] 

I 1/129 

sin^ 9w 

0.230 

Gf 

1.166 X lO-^GeV^ 

Vud 

0.97427 

Vus 

0.22534 

Re(r) 

0.001513 

Im(T) 

-0.000601 


TABLE IV: 'Zi and y^. The parameters of the standard model tabulated in Table Hill 

are used in the calculations. We set y = uric = 1.3 GeV, and choose two values q* = 1/a and vr/a 
as the matching scale from the lattice to the continuum. 













q* 

= 

1/a 






q* 

= 

tt/o 




i 










Zi 




Vi 




Zi 




Vi 

1 - 

-4.184 

X 

10- 

-1 




0 

-4.487 

X 

10- 

-1 




0 - 

-3.453 

X 

10- 

-1 




0 

2 

1.218 

X 

lO+o 




0 

1.029 

X 

lO+o 




0 

9.321 

X 

10- 

-1 




0 

3 

4.575 

X 

10- 

-3 

2.910 

X 

10- 

-2 

-9.327 

X 

10- 

-3 

3.145 

X 

10' 

-2 

-1.246 

X 

10- 

-2 

2.629 

X 

10' 

-2 

4 - 

-1.373 

X 

10- 

-2 

-5.782 

X 

10- 

-2 

2.703 

X 

10- 

-2 

-5.628 

X 

10' 

-2 

3.173 

X 

10- 

-2 

-5.108 

X 

10' 

-2 

5 

4.575 

X 

10- 

-3 

4.869 

X 

10- 

-3 

-7.309 

X 

10- 

-3 

1.402 

X 

10' 

-2 

-9.284 

X 

10- 

-3 

1.470 

X 

10' 

-2 

6 - 

-1.373 

X 

10- 

-2 

-9.009 

X 

10- 

-2 

2.323 

X 

10- 

-2 

-4.700 

X 

10' 

-2 

2.130 

X 

10- 

-2 

-4.015 

X 

10' 

-2 

7 

6.305 

X 

10- 

-5 

-2.010 

X 

10- 

-4 

5.777 

X 

10- 

-5 

-2.514 

X 

10' 

-4 

1.731 

X 

10- 

-4 

-1.904 

X 

10' 

-4 

8 




0 

1.098 

X 

10- 

-3 

-4.573 

X 

10- 

-5 

4.600 

X 

10' 

-4 . 

-6.871 

X 

10- 

-5 

3.183 

X 

10' 

-4 

9 

6.305 

X 

10- 

-5 

-1.168 

X 

10- 

-2 

-3.047 

X 

10- 

-6 

-9.925 

X 

10' 

-3 

7.287 

X 

10- 

-5 

-8.995 

X 

10' 

-3 

10 




0 

4.357 

X 

10- 

-3 

5.277 

X 

10- 

-5 

4.635 

X 

10' 

-3 

6.368 

X 

10- 

-5 

3.717 

X 

10' 

-3 
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TABLE V: Results of the K —)• tttt decw amplitudes. The results by the RBC-UKQCD Collab¬ 
oration at m-jr = 422 MeV 3|, 330MeV (^, the physical quark mass in the continuum limit (only 
for the A/ = 3/2 process) (3|, and the experimental values are also tabulated. 



q* = 1/a 

q* = tt/ a 


RBG-UKQGD 

Exp. 

a (fm) 

0.091 


0.114 

0.114 

- 

TTljr (MeV) 

280 


330 

422 

140 140 

ReA 2 (xlO-® GeV) 

2.426(38) 

2.460(38) 

2.668(14) 

4.911(31) 

1.50(4)(14) 1.479(4) 

ReAo(xlO-® GeV) 

60(36) 

56(32) 

31.1(4.5) 

38.0(8.2) 

33.2(2) 

ReAo /ReA2 

25(15) 

23(13) 

12.0(1.7) 

7.7(1.7) 

22.45(6) 

ImA2(x 10-12 Qgv) 

-1.192(14) -0.7457(83) 

-0.6509(34) 

-0.5502(40) - 

0.699(20) (84) 

ImAo(x 10-12 GeV) 

-67(56) 

-52(48) 

-33(15) 

-25(22) 


Re(eVe)(xlO-3) 

0.8(2.5) 

0.9(2.5) 

2.0(1.7) 

2.7(2.6) 

1.66(23) 
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